* Begin log file 
capture log close sublog
log using "$projdir/log/8_figures.txt", name(sublog) text replace

*-------------------------------------------------------------------------------
* Description: Make the figures for the main text
* Authors: Justin Wiltshire, Carl McPherson, Denis Sosinksiy
*
* Updated: May 17, 2024
*-------------------------------------------------------------------------------
clear

* Generate tempfiles
tempfile core
tempfile core0
tempfile core1
tempfile core2

* Ensure necessary folders exist
capture mkdir "$projdir/out/figures"

*-------------------------------------------------------------------------------
* Figure 1:
* Distribution of Average Wage by County
*-------------------------------------------------------------------------------

* Set start and end quarters
local start_tm = 199
local end_tm = 239
local tot_tm = `end_tm' - `start_tm' + 1
local endqtrlab "2019Q4"
if `end_tm' == 248 {
	local endqtrlab "2022Q1"
}

* Ensure distinct package is installed
capture ssc install distinct

* Load the data
use "$projdir/dta/build/cln/analysis_panel_cty.dta", clear

keep if !mi(countyfips) & !mi(avg_wkly_wage) & !mi(employment) & avg_wkly_wage != 0 & employment != 0
summ employment if ind == 10 & year == 2015, d 

collapse (rawsum) employment (min) min_emp (mean) avg_wkly_wage [fw=round(employment)], by(year statef county* area* ind)

summ avg_wkly_wage [fw=round(employment)] if year == 2013 & ind == 10, d 
local p10 = r(p10)
local p25 = r(p25)
local mean = r(mean)
local p75 = r(p75)
local p90 = r(p90)

* Generate donor variable
gen donor = 1
replace donor = 0  if statef == 6 
replace donor = 0  if statef == 36 //NY

gen mid = inlist(area,36085, 36087, 36027, 36091, 36065, 36063, 36007, 36111, 36083, 36093, 36069)
replace donor = 0 if mid == 1 
keep if donor == 0

drop if inlist(area, 36029, 6079)

drop if mid == 1

qui gen x = (ind == 722 & min_emp >= 5000)
bysort statefips countyfips: egen empkeep = max(x)
keep if empkeep == 1

gen name = subinstr(subinstr(subinstr(area_title, " County", "",.), "New York","NY",.), "California", "CA", .)  
replace name = "New York, NY" if area == 36061

drop if !inlist(statef,6,36)

keep if year == 2013
keep if ind == 10 

gen localmw = inlist(area, 6001, 6013, 6037, 6041, 6073, 6075, 6081, 6085, 6097) ///
				| inlist(area, 36005,36047,36061,36081,36085)				

gen ny_large_no = avg_wkly_wage if statef == 36 & mid == 0 & localmw == 0
gen ny_large_yes = avg_wkly_wage if statef == 36 & mid == 0 & localmw == 1
gen ca_noloc = avg_wkly_wage if statef == 6 & localmw == 0
gen ca_yesloc = avg_wkly_wage if statef == 6 & localmw == 1

foreach v in ca_yesloc ca_noloc ny_large_no ny_large_yes {
	di "`v'"
	tab `v'
	di ""
}

#delimit ;
graph dot ca_yesloc ca_noloc ny_large_yes ny_large_no
	, 
	over(name, sort(avg_wkly_wage) lab(labsize(vsmall))) 
	linetype(line) lines(lw(vvthin) lp(blank))
	
	marker(1, msize(small) mcolor(navy)) 
	marker(2, msize(small) ms(Oh) mcolor(navy)) 
	marker(3, msize(small) ms(T) mcolor(cranberry)) 
	marker(4, msize(small) ms(Th) mcolor(cranberry)) 
	
	yline(`p10', lcolor(red) lp(dash))
	yline(`p25', lcolor(red) lp(solid)) 
	yline(`mean', lcolor(black) lp(solid) lwidth(medthick))
	yline(`p75', lcolor(red) lp(solid)) 
	yline(`p90', lcolor(red) lp(dash))
	
	legend(order(1 "California, Local" 
				 2 "California, No Local"
				 3 "New York, Local"
				 4 "New York, No Local"
				 ) pos(2) ring(0) row(6) size(small))

	ytitle("Average weekly wages") ylabel(0(500)2250) ymtick(##2)
	graphregion(color(white)) bgcolor(white)
	outergap(0)
	
	xsize(5.4)
	ysize(6.6)
;
#delimit cr 

* Save the graph
graph export "$projdir/out/figures/figure1.pdf", replace 

*-------------------------------------------------------------------------------
* Figure 2 (upper panel):
* Bites for California
*-------------------------------------------------------------------------------
use "$projdir/dta/build/cln/epi_cps_panel", clear 
keep if statef == 6  

* Rename wage hrwage
rename orgw earnwt

* Get upcoming mw
egen group = group(max_mw)

preserve

	collapse (first) max_mw (last) quarterly_date, by(group)
	sort quarterly_date
	gen next_mw = max_mw[_n+1] 
	replace next_mw = 15 if next_mw == . 
	
	keep group next_mw
	
	tempfile next 
	save `next'
	
restore 

merge m:1 group using `next', assert(3) nogen
replace next_mw = 15.50 if year == 2022

* Keep only people in the labor force and non-self-emploted
drop if lfstat != 1
drop if selfemp != 0 

* Fraction under next 
gen under = (hrwage < next_mw)*100
replace under = . if hrwage == . 

replace next_mw = 10 if year == 2014 
replace max_mw = 9 if year == 2014 

preserve 

* Collapse to qtr-year 
#delimit ;
collapse 
	(median) p50_hw = hrwage  
	(mean) avg_hw = hrwage
		   under
	[w=earnwt]
	,
	by(statefip year max_mw next_mw)
;
#delimit cr 

gen countyfips = -1

tempfile bstate 
save `bstate'

restore 

gen cnt =  1 

preserve 

* Collapse to county 
#delimit ;
collapse 
	(median) p50_hw = hrwage  
	(mean) avg_hw = hrwage
		   under
	[w=earnwt]
	,
	by(statefip countyfip year max_mw next_mw cnt)
;
#delimit cr

tempfile bcounty
save `bcounty'

restore 

* Food 
gen  food = inlist(ind17,8680,8690) | inlist(ind12,8680,8690) | inlist(ind07,8680,8690)
forv n = 15/23 {
gen under`n' = hrwage < `n'
}
drop if hrwage == . 

format quarterly_date %tq 

*Collapse to qtr-year 
#delimit ;
collapse 
	(median) p50_hw = hrwage  
	(mean) avg_hw = hrwage
		   under
	[w=earnwt] if food == 1 
	,
	by(statefip year max_mw next_mw)
;
#delimit cr 

gen countyfips = -3

append using `bstate'
append using `bcounty'

*Get ratio 
gen kratio = max_mw / p50_hw*100 

global xline "xline(2013.5, lcolor(gs10) lp(dash))"

local ystart = 2010
local yend = 2022
keep if inrange(year, `ystart', `yend')
sort year 

*br year kratio under if countyfips == -1 

* Plot ratio 
#delimit ;
tw 	(scatter kratio year if countyfips == -1, c(l))
	(scatter kratio year if countyfips == 75, c(l) lcolor(cranberry) mcolor(cranberry) lp(longdash) ms(S))
	(scatter kratio year if countyfips == 19, c(l) lcolor(green) mcolor(green) lp(dash) ms(T))
	(scatter kratio year if countyfips == -3, c(l) lcolor(gold) mcolor(gold) lp(dash_dot) ms(D))
	,
	xtitle("") xlabel(`ystart'(1)`yend', angle(45)) xtick(`ystart'(1)`yend')
	ytitle("Minimum / median hourly wage")
	ylabel(0(10)100)
	legend(lab(1 "All California") 
		   lab(2 "San Francisco") 
		   lab(3 "Fresno")
		   lab(4 "Restaurants")
		   row(2) pos(6) span)
	${xline}
	graphregion(color(white)) bgcolor(white)

;
#delimit cr  
graph export "$projdir/out/figures/figure2a.pdf", replace 

*-------------------------------------------------------------------------------
* Figure 2 (lower panel):
* Fraction earning under upcoming MW
*-------------------------------------------------------------------------------
#delimit ;
tw 	(scatter under year if countyfips == -1, c(l))
	(scatter under year if countyfips == 75, c(l) lcolor(cranberry) mcolor(cranberry) lp(longdash) ms(S))
	(scatter under year if countyfips == 19, c(l) lcolor(green) mcolor(green) lp(dash) ms(T))
	(scatter under year if countyfips == -3, c(l) lcolor(gold) mcolor(gold) lp(dash_dot) ms(D))
	,
	xtitle("") xlabel(`ystart'(1)`yend', angle(45)) xtick(`ystart'(1)`yend')
	ytitle("% workers paid under upcoming min wage")
	ylabel(0(10)60)
	legend(lab(1 "All California") 
		   lab(2 "San Francisco") 
		   lab(3 "Fresno")
		   lab(4 "Restaurants")
		   row(2) pos(6) span)
	${xline}
	graphregion(color(white)) bgcolor(white)

;
#delimit cr  

graph export "$projdir/out/figures/figure2b.pdf", replace 

*-------------------------------------------------------------------------------
* Figure 3 (upper panel):
* Average and marginal effects figures (all counties, Bias-corrected)
*-------------------------------------------------------------------------------

* Make the figures
foreach y in avg_wkly_wage employment {
	local n = 722513
	
	* Get the ATT
	qui use "$projdir/dta/analysis/qcew/`y'/naics`n'/agg_att/bias_corrected/bias_corr_att.dta", clear
	qui rename _tm evt_time
	qui gen cty_fips = 0
	qui save "`core'", replace
	
	* Get the 2010 population weights for each treated county
	foreach z in California NY {
		
		* Get filelist
		local filelist: dir "$projdir/dta/analysis/qcew/`y'/naics`n'/`z'/corrected" files "*.dta"
		if "`z'" == "California" {
			qui gettoken first filelist : filelist
			qui use "$projdir/dta/analysis/qcew/`y'/naics`n'/`z'/corrected/`first'", clear
		}

		* Append the rest of the data
		foreach c in `filelist' {
			if "`c'" != "`first'" {
				qui append using "$projdir/dta/analysis/qcew/`y'/naics`n'/`z'/corrected/`c'"
			}
		}
	}
	
	* Clean and save tempfile
	qui keep if cty_fips == trunit
	bysort cty_fips: keep if _n == 1
	qui keep cty_fips pop10
	qui compress
	qui save "`core2'", replace
	
	* Append the bias-corrected data files
	foreach z in California NY  {
		
	* Get filelist
	local filelist: dir "$projdir/dta/analysis/qcew/`y'/naics`n'/`z'" files "*.dta"
	if "`z'" == "California" {
			qui gettoken first filelist : filelist
			qui use "$projdir/dta/analysis/qcew/`y'/naics`n'/`z'/`first'", clear
		}
		
	* Append the rest of the data
	foreach c in `filelist' {
		if "`c'" != "`first'" {
			qui append using "$projdir/dta/analysis/qcew/`y'/naics`n'/`z'/`c'"
			}
		}
	}
	
	* Clean data
	qui keep if cty_fips == trunit
	qui drop if mi(_Y_treated)
	
	* Merge in the 2010 population weights and keep key variables
	merge m:1 cty_fips using "`core2'", nogen norep keep(3)
	
	* Exclude mid counties in NY
	drop if inlist(cty_fips, 36007, 36027, 36063, 36065, 36069, 36083, 36085, 36087, 36091, 36093, 36111) 
		
	* Put in event time and keep a sample balanced in event time
	qui gen evt_time = _time - trperiod
	qui keep if inrange(evt_time, -17, 21)
			
	* Keep key variables
	qui keep cty_fips gap_bc evt_time pop10
	
	* Save the tempfile
	qui compress
	save "`core2'", replace

	* Append the ATT
	qui append using "`core'"
	qui keep if inrange(evt_time, -17, 21)
	
	* Exclude mid counties in NY
	drop if inlist(cty_fips, 36007, 36027, 36063, 36065, 36069, 36083, 36085, 36087, 36091, 36093, 36111) 
		
	* Set macros
	local ytitle "Difference (percent)"
	local tit0
	local ymin = -40
	local ymax = 45
	local yint = 10
	if strpos("`y'", "avg_wkly_wage") != 0 {
		local tit "Average weekly earnings"
		local nam "l"
	}
	if strpos("`y'", "employment") != 0 {
		local tit "Employment"
		local nam "r"
	} 	
	
	tsset cty_fips evt_time
	#delimit ;
		twoway  
			(scatter gap_bc evt_time [aw=pop10], mcolor(eltblue) msymbol(circle_hollow) mlwidth(vvthin) msize(tiny))			
			(line gap_bc evt_time if cty_fips == 0, lwidth(medthick) lcolor(blue) lp(solid)),
			xline(0, lp(dot) lcolor(gs8))
			yline(0, lp(longdash) lcolor(black))
			subtitle("`tit0'`tit'")
			xtitle("Event Quarter") ytitle("`ytitle'")
			yscale(range(`ymin' `ymax'))
			ylabel(`ymin'(`yint')`ymax')
			xlabel(-16(2)20, angle(45))
			xscale(range(-17 21))
			plotregion(margin(zero))
			graphregion(fcolor(white))
			legend(rows(2) order(2 "Average county-level treatment effect" 1 "Individual county treatment effects"));
		#delimit cr
		
	* Save the graph
	graph export "$projdir/out/figures/figure3a`nam'.pdf", replace
}

*-------------------------------------------------------------------------------
* Figure 3 (lower panel): 
* Average county figures (all counties, Bias-corrected)
*-------------------------------------------------------------------------------

* Make the figures
foreach y in employment avg_wkly_wage {
	local n = 722513
	
	* Load the data
	qui use "$projdir/dta/analysis/qcew/`y'/naics`n'/agg_att/bias_corrected/bias_corr_att_distn.dta", clear
	
	* Keep time periods of interest
	qui keep if inrange(_tm, -17, 21)
	
	* Set macros
	local tit
	local ytitle "Difference (percent)"
	local tit0
	local ymin = -20
	local ymax = 25
	local yint = 5
	if strpos("`y'", "avg_wkly_wage") != 0 {
			local tit "Average weekly earnings"
			local nam "l"
	}
	if strpos("`y'", "employment") != 0 {
		local tit "Employment"
		local nam "r"
	} 
	
	* Set the plot lists
	qui distinct _placeboID
	local ciid = r(ndistinct) + 1
	local trid = r(ndistinct) + 2
	qui levelsof _placeboID, local(id)
	local plotlist
	
	foreach c of local id {
		local plotlist "`plotlist' (line gap_bc _tm if _placeboID == `c' & inrange(gap_bc, `ymin', `ymax'), lwidth(vthin) lcolor(gs8) lp(solid))"
	}
	
				
	di "Graph the distribution of estimated classic treatment effect paths for `y'"	
	tsset _placeboID _tm
	#delimit ;
		twoway  
			`plotlist'
			(rarea LB_95 UB_95 _tm if _placeboID == 0, color(gs12%50) lwidth(none))
			(line gap_bc _tm if _placeboID == 0, lwidth(medthick) lcolor(blue) lp(solid)),
			xline(0, lp(dot) lcolor(gs8))
			yline(0, lp(longdash) lcolor(black))
			subtitle("`tit0'`tit'")
			xtitle("Event Quarter") ytitle("`ytitle'")
			yscale(range(`ymin' `ymax'))
			ylabel(`ymin'(`yint')`ymax')
			xlabel(-16(2)20, angle(45))
			xscale(range(-17 21))
			plotregion(margin(zero))
			graphregion(fcolor(white))
			legend(order(`trid' "Treated county ATT" 1 "Sample placebo ATTs" `ciid' "95% Confidence intervals"));
		#delimit cr

	* Save the graphs
	graph export "$projdir/out/figures/figure3b`nam'.pdf", replace
}

*-------------------------------------------------------------------------------
* Figure 4 (upper panel): 
* Relative employment level for treated vs donor counties
*-------------------------------------------------------------------------------	

* Set start and end quarters
local trp = 218
local start_tm = 199
local end_tm = 251 
local tot_tm = `end_tm' - `start_tm' + 1
local endqtrlab "2019Q4"
if `end_tm' == 251 {
	local endqtrlab "2022Q4"
}

* Ensure distinct package is installed
capture ssc install distinct

* Load the data
use "$projdir/dta/build/cln/analysis_panel_cty.dta", clear

* Rename key variables
qui rename (quarterly_date earnings lnemp ind) (tm demean_lnearn demean_lnemp naics)
format tm %tq

* Generate unique county id and drop obs which don't observe it
qui gen cty_fips = 1000*statefips + countyfips
qui keep if !mi(cty_fips) & !mi(avg_wkly_wage) & !mi(employment) & avg_wkly_wage != 0 & employment != 0
assert cty_fips
assert avg_wkly_wage
assert employment
assert !mi(avg_wkly_wage)
assert !mi(employment)

* Recode pre-2011 full- and limited-service restaurant NAICS codes
qui recode naics (722211=722513) if year <= 2011

* Keep NAICS 10, 722, and 722513
qui keep if inlist(naics, 10, 722, 722513)

* Sort
sort cty_fips tm naics

* Restrict to desired time-period
keep if inrange(tm, `start_tm', `end_tm')

* Identify treated and donor counties
levelsof cty_fips if inlist(statefips, 6, 36), local(trcty)
local trctylist = subinstr("`trcty'", " ", ",",.)
qui distinct cty_fips if inlist(cty_fips, `trctylist')
local trcty_count = r(ndistinct)

* Generate donor variable
qui gen donor = !inlist(cty_fips, `trctylist')

* Restrict to donor counties with no min wage changes (and the treated counties)
qui egen double max_loc_max = max(loc_max_mw), by(cty_fips)
qui egen double min_loc_max = min(loc_max_mw), by(cty_fips)
qui gen x = (naics == 722 & min_emp >= 5000)
qui drop if cty_fips == 13121
bysort cty_fips: egen empkeep = max(x)
qui keep if (max_loc_max == min_loc_max | donor == 0) & empkeep == 1
if `end_tm' > 239 {
	qui drop if inlist(floor(cty_fips/1000), 25, 51)
}

* Save tempfile
qui compress
qui save "`core'", replace

* Make county change adjustments and get mean of avg_wkly_wage over remaining counties
qui do "$projdir/do/build/cw/cty_changes_NYCnocombine.do"
qui collapse (mean) avg_wkly_wage donor year qtr [aw=employment], by(cty_fips tm naics)
assert inlist(donor, 0, 1)

* Save tempfile
qui save "`core2'", replace

* Reload the data
qui use "`core'", clear

* Make county change adjustments and get sum of employment over remaining counties
qui do "$projdir/do/build/cw/cty_changes_NYCnocombine.do"
qui collapse (sum) employment pop10, by(cty_fips tm naics)

* Merge
qui merge 1:1 cty_fips tm naics using "`core2'", nogen norep

* Keep only a balanced panel over periods of interest
qui keep if inrange(tm, `start_tm', `end_tm')	
bysort cty_fips naics: gen ct = _N
qui keep if ct == `tot_tm'

* Merge with the LAUS data
qui merge m:1 cty_fips year using "$projdir/dta/build/cln/laus_cty_NYCnocombine.dta", nogen norep keep(3)

* Merge with the Covid-impact index
qui save "`core'", replace
qui use "$projdir/dta/build/cln/covid_index.dta", clear
qui rename (countyfips) (cty_fips)
qui gen covid_index = (workplace + retail)/2
qui keep cty_fips covid_index
qui save "`core2'", replace
qui use "`core'", clear
merge m:1 cty_fips using "`core2'", nogen norep keep(3)

* Observe overall employment in each period
qui gen x = employment if naics == 10
bysort cty_fips tm: egen emp10 = max(x)
qui drop x
qui drop if naics == 10

* Keep variables of interest	
qui keep cty_fips tm naics avg_wkly_wage employment emp10 pop10 unr covid_index donor

* Save
qui save "`core0'", replace

* For each of 722 and 722513
foreach n in 722 722513 {
	
	* Load and restrict
	qui use "`core0'", clear
	qui keep if naics == `n'

	* Put in % of 2014q2 levels
	foreach v in employment avg_wkly_wage emp10 {
		qui gen x = `v' if tm == 217
		bysort cty_fips: egen xx = max(x)
		qui replace `v' = 100*`v'/xx
		qui drop x xx
	}

	* Reshape data to debias it
	local predvars ""
	forval t = 199/207 {
		qui gen x = employment if tm == `t'
		bysort cty_fips: egen emp_`t' = max(x)
		qui drop x
		local predvars "`predvars' emp_`t'"
		qui gen x = avg_wkly_wage if tm == `t'
		bysort cty_fips: egen earn_`t' = max(x)
		qui drop x
		local predvars "`predvars' earn_`t'"
		qui gen x = emp10 if tm == `t'
		bysort cty_fips: egen emp10_`t' = max(x)
		qui drop x
		local predvars "`predvars' emp10_`t'"
	}
	foreach v in employment avg_wkly_wage emp10 unr {
		bysort cty_fips: egen x = mean(`v') if inrange(tm, 199, 207)
		bysort cty_fips: egen mean_`v' = max(x)
		qui drop x
		local predvars "`predvars' mean_`v'"
	}
	save "`core1'", replace

	* "Covid- and bias-correct" the outcomes of interest
	qui gen p_employment = .
	qui gen p_avg_wkly_wage = .
	qui levelsof tm, local(Time)
	foreach v in employment avg_wkly_wage {
		foreach t of local Time {
			qui reg `v' `predvars' covid_index [aw=pop10] if tm == `t' & donor == 1
			qui predict p_`v'_`t', resid
			qui replace p_`v' = p_`v'_`t' if tm == `t'
			qui drop p_`v'_`t'
		}
	}
		
	* Keep key variables
	qui keep cty_fips tm naics employment avg_wkly_wage p_* emp10 pop10 unr donor
		
	* Compress, sort, and save
	qui compress
	qui sort cty_fips tm
	if `n' == 722513 {
		qui append using "`core'"
	}
	qui save "`core'", replace
}

* Plot the donor pool Covid-corrected earnings
qui use "`core'", clear
qui keep if naics == 722513

* Collapse by period
collapse (mean) avg_wkly_wage employment [aw=pop10], by(tm donor)
		
* Append to tempfile
gen cty_fips = 0
qui append using "`core'"
qui keep if naics == 722513 | cty_fips == 0
qui distinct cty_fips if donor == 1
local treated = r(ndistinct) + 1
qui distinct cty_fips if cty_fips != 0
local donoravg = r(ndistinct) + 1
local travg = r(ndistinct) + 2

local y = "employment" 

* Set the plot lists for the ATT and the full distribution of estimated county-level treatment effects
qui levelsof cty_fips if donor == 1 & cty_fips != 0, local(id)
local donorplotlist
local trplotlist
foreach c of local id {
	local donorplotlist "`donorplotlist' (line `y' tm if cty_fips == `c', lcolor(gs8) lwidth(vthin) lp(solid))"
}
qui levelsof cty_fips if donor == 0 & cty_fips != 0, local(id)
foreach c of local id {
	local trplotlist "`trplotlist' (line `y' tm if cty_fips == `c', lcolor(eltblue) lwidth(vthin) lp(solid))"
}
		
* Set macros
local ytitle "Relative employment level (percent)"
local tit0
local ymin = 70
local ymax = 130
local yint = 10

* Figure 3 Panel A
#delimit ;
	twoway  
		(line `y' tm if cty_fips == 0 & donor == 0, lcolor(blue) lwidth(medthick) lp(solid))
		(line `y' tm if cty_fips == 0 & donor == 1, lcolor(black) lwidth(medium) lp(dash))
		,
		xline(218, lp(dot) lcolor(gs8))
		yline(100, lp(longdash) lcolor(black))
		xtitle(" ") ytitle("`ytitle'")
		yscale(range(`ymin' `ymax'))
		ylabel(`ymin'(`yint')`ymax')
		tlabel(2009Q4(4)2022Q4, angle(45))
		plotregion(margin(zero))
		graphregion(fcolor(white))
		legend(rows(4) order(1 "Treated county average" 2 "Donor county average"));
	#delimit cr
				
* Save the graph
graph export "$projdir/out/figures/figure4a.pdf", replace	

*-------------------------------------------------------------------------------
* Figure 4 (lower left panel):
* Average county figure (all countuies, Bias-corrected, through e33)
*-------------------------------------------------------------------------------

* Make the figures
local y  = "employment"
local n = 722513
	
* Load the data
qui use "$projdir/dta/analysis/qcew/`y'/naics`n'/agg_att/bias_corrected/bias_corr_att_distn.dta", clear
	
* Set macros
local tit
local ytitle "Difference (percent)"
local tit0
local ymin = -20
local ymax = 20
local yint = 5
if strpos("`y'", "employment") != 0 {
	local tit "Uncorrected"
} 
	
* Set the plot lists
qui distinct _placeboID
local ciid = r(ndistinct) + 1
local trid = r(ndistinct) + 2
qui levelsof _placeboID, local(id)
local plotlist
	
foreach c of local id {
	local plotlist "`plotlist' (line gap_bc _tm if _placeboID == `c' & inrange(gap_bc, `ymin', `ymax'), lwidth(vthin) lcolor(gs8) lp(solid))"
}
				
tsset _placeboID _tm
#delimit ;
	twoway  
		`plotlist'
		(rarea LB_95 UB_95 _tm if _placeboID == 0, color(gs12%50) lwidth(none))
		(line gap_bc _tm if _placeboID == 0, lwidth(medthick) lcolor(blue) lp(solid)),
		xline(0, lp(dot) lcolor(gs8))
		yline(0, lp(longdash) lcolor(black))
		subtitle("`tit0'`tit'")
		xtitle("Event Quarter") ytitle("`ytitle'")
		yscale(range(`ymin' `ymax'))
		ylabel(`ymin'(`yint')`ymax')
		xlabel(-16(2)32, angle(45))
		xscale(range(-17 33))
		plotregion(margin(zero))
		graphregion(fcolor(white))
		legend(order(`trid' "Treated county ATT" 1 "Sample placebo ATTs" `ciid' "95% Confidence intervals"));
	#delimit cr

* Save the graph
graph export "$projdir/out/figures/figure4bl.pdf", replace

*-------------------------------------------------------------------------------
* Figure 4 (lower right panel): 
* Average county figure (all counties, Covid-corrected, through e33)
*-------------------------------------------------------------------------------

* Set locals
local y = "employment"
local n = 722513
	
* Load the data
qui use "$projdir/dta/analysis/qcew/`y'/naics`n'/agg_att/covid_corrected/covid_corr_att_distn.dta", clear
	
* Keep time periods of interest
qui keep if inrange(_tm, -17, 33)
	
* Set macros
local tit
local ytitle "Difference (percent)"
local tit0
local ymin = -20
local ymax = 20
local yint = 5
if strpos("`y'", "avg_wkly_wage") != 0 {
		local tit "Average weekly earnings"
}
if strpos("`y'", "employment") != 0 {
	local tit "Covid-corrected"
} 
	
* Set the plot lists
qui distinct _placeboID
local ciid = r(ndistinct) + 1
local trid = r(ndistinct) + 2
qui levelsof _placeboID, local(id)
local plotlist
	
foreach c of local id {
	local plotlist "`plotlist' (line p_gap _tm if _placeboID == `c' & inrange(p_gap, `ymin', `ymax'), lwidth(vthin) lcolor(gs8) lp(solid))"
}
				
di "Graph the distribution of estimated classic treatment effect paths for `y'"	
tsset _placeboID _tm
#delimit ;
	twoway  
		`plotlist'
		(rarea LB_95 UB_95 _tm if _placeboID == 0, color(gs12%50) lwidth(none))
		(line p_gap _tm if _placeboID == 0, lwidth(medthick) lcolor(blue) lp(solid)),
		xline(0, lp(dot) lcolor(gs8))
		yline(0, lp(longdash) lcolor(black))
		subtitle("`tit0'`tit'")
		xtitle("Event Quarter") ytitle("`ytitle'")
		yscale(range(`ymin' `ymax'))
		ylabel(`ymin'(`yint')`ymax')
		xlabel(-16(2)32, angle(45))
		xscale(range(-17 33))
		plotregion(margin(zero))
		graphregion(fcolor(white))
		legend(order(`trid' "Treated county ATT" 1 "Sample placebo ATTs" `ciid' "95% Confidence intervals"));
	#delimit cr

* Save the graph
graph export "$projdir/out/figures/figure4br.pdf", replace


*-------------------------------------------------------------------------------
* Figure 5 (top panel): 
* Change In Time Spent in Restaurants and Retail Establishments
*-------------------------------------------------------------------------------
* Load Google data
use "$projdir/dta/build/src/gm_state.dta", clear 

gen date = mdy(month, day, year)
format %td date

gen w = week(date)
gen weekly = yw(year, w)
format %tw weekly

bysort statefips weekly: egen last_day_w = max(date)
format %tw last_day_w
by statefips weekly: keep if date == last_day_w

rename gps_retail_and_recreation retail
rename gps_workplaces workplace

keep statefips retail workplace date weekly

sort statefips weekly

generate donor_st = 0
replace donor_st = 1 if ///
inlist(statefips, 1,13,16,18,19,20,21,22,28,33,37,38,40,42,45,47,48,49,55,56)

cap drop upper 
gen upper = 0.3

#delimit ;
tw 	
	(area upper date if inrange(date, 21989,22111), fcolor(gs14*0.6) lcolor(gs14*0.6) base(-.7))
	(scatter retail date  if statefips == 1, c(l) lcolor(gs10) mcolor(gs10) ms(none) lp(solid))
	(scatter retail date  if statefips == 13, c(l) lcolor(gs10) mcolor(gs10) ms(none) lp(solid))
	(scatter retail date  if statefips == 16, c(l) lcolor(gs10) mcolor(gs10) ms(none) lp(solid))
	(scatter retail date  if statefips == 18, c(l) lcolor(gs10) mcolor(gs10) ms(none) lp(solid))
	(scatter retail date  if statefips == 19, c(l) lcolor(gs10) mcolor(gs10) ms(none) lp(solid))
	(scatter retail date  if statefips == 20, c(l) lcolor(gs10) mcolor(gs10) ms(none) lp(solid))
	(scatter retail date  if statefips == 21, c(l) lcolor(gs10) mcolor(gs10) ms(none) lp(solid))
	(scatter retail date  if statefips == 22, c(l) lcolor(gs10) mcolor(gs10) ms(none) lp(solid))
	(scatter retail date  if statefips == 28, c(l) lcolor(gs10) mcolor(gs10) ms(none) lp(solid))
	(scatter retail date  if statefips == 33, c(l) lcolor(gs10) mcolor(gs10) ms(none) lp(solid))
	(scatter retail date  if statefips == 37, c(l) lcolor(gs10) mcolor(gs10) ms(none) lp(solid))
	(scatter retail date  if statefips == 38, c(l) lcolor(gs10) mcolor(gs10) ms(none) lp(solid))
	(scatter retail date  if statefips == 40, c(l) lcolor(gs10) mcolor(gs10) ms(none) lp(solid))
	(scatter retail date  if statefips == 42, c(l) lcolor(gs10) mcolor(gs10) ms(none) lp(solid))
	(scatter retail date  if statefips == 45, c(l) lcolor(gs10) mcolor(gs10) ms(none) lp(solid))
	(scatter retail date  if statefips == 47, c(l) lcolor(gs10) mcolor(gs10) ms(none) lp(solid))
	(scatter retail date  if statefips == 48, c(l) lcolor(gs10) mcolor(gs10) ms(none) lp(solid))
	(scatter retail date if statefips == 49, c(l) lcolor(gs10) mcolor(gs10) ms(none) lp(solid))
	(scatter retail date if statefips == 55, c(l) lcolor(gs10) mcolor(gs10) ms(none) lp(solid))
	(scatter retail date if statefips == 56, c(l) lcolor(gs10) mcolor(gs10) ms(none) lp(solid))
	(scatter retail date  if statefips == 6, c(l) lcolor(blue) ms(none) lwidth(medthick) lp(solid))
	(scatter retail date  if statefips == 36, c(l) lcolor(cranberry) ms(none) lwidth(medthick) lp(dash))
	,
	xtitle("") xmtick(##2)
	ytitle("Percent change in time spent in restaurants/retail")
	legend(order(22 "California" 23 "New York" 2 "Donor States") row(1) span) 
	ylab(-0.7(0.2)0.3)
	yline(0, lcolor(black))
	graphregion(color(white)) bgcolor(white)
;
#delimit cr 
graph export "$projdir/out/figures/figure5a.pdf",  replace

*-------------------------------------------------------------------------------
* Figure 5 (bottom panel): 
* Change In Time Spent in Workplace Location
*-------------------------------------------------------------------------------
cap drop upper 
gen upper = 0.1
#delimit ;
tw 	
	(area upper date if inrange(date, 21989,22111), fcolor(gs14*0.6) lcolor(gs14*0.6) base(-.6))
	(scatter workplace date  if statefips == 1, c(l) lcolor(gs10) mcolor(gs10) ms(none) lp(solid))
	(scatter workplace date  if statefips == 13, c(l) lcolor(gs10) mcolor(gs10) ms(none) lp(solid))
	(scatter workplace date  if statefips == 16, c(l) lcolor(gs10) mcolor(gs10) ms(none) lp(solid))
	(scatter workplace date  if statefips == 18, c(l) lcolor(gs10) mcolor(gs10) ms(none) lp(solid))
	(scatter workplace date  if statefips == 19, c(l) lcolor(gs10) mcolor(gs10) ms(none) lp(solid))
	(scatter workplace date  if statefips == 20, c(l) lcolor(gs10) mcolor(gs10) ms(none) lp(solid))
	(scatter workplace date  if statefips == 21, c(l) lcolor(gs10) mcolor(gs10) ms(none) lp(solid))
	(scatter workplace date  if statefips == 22, c(l) lcolor(gs10) mcolor(gs10) ms(none) lp(solid))
	(scatter workplace date  if statefips == 28, c(l) lcolor(gs10) mcolor(gs10) ms(none) lp(solid))
	(scatter workplace date  if statefips == 33, c(l) lcolor(gs10) mcolor(gs10) ms(none) lp(solid))
	(scatter workplace date  if statefips == 37, c(l) lcolor(gs10) mcolor(gs10) ms(none) lp(solid))
	(scatter workplace date  if statefips == 38, c(l) lcolor(gs10) mcolor(gs10) ms(none) lp(solid))
	(scatter workplace date  if statefips == 40, c(l) lcolor(gs10) mcolor(gs10) ms(none) lp(solid))
	(scatter workplace date  if statefips == 42, c(l) lcolor(gs10) mcolor(gs10) ms(none) lp(solid))
	(scatter workplace date  if statefips == 45, c(l) lcolor(gs10) mcolor(gs10) ms(none) lp(solid))
	(scatter workplace date  if statefips == 47, c(l) lcolor(gs10) mcolor(gs10) ms(none) lp(solid))
	(scatter workplace date  if statefips == 48, c(l) lcolor(gs10) mcolor(gs10) ms(none) lp(solid))
	(scatter workplace date if statefips == 49, c(l) lcolor(gs10) mcolor(gs10) ms(none) lp(solid))
	(scatter workplace date if statefips == 55, c(l) lcolor(gs10) mcolor(gs10) ms(none) lp(solid))
	(scatter workplace date if statefips == 56, c(l) lcolor(gs10) mcolor(gs10) ms(none) lp(solid))
	(scatter workplace date  if statefips == 6, c(l) lcolor(blue) ms(none) lwidth(medthick) lp(solid))
	(scatter workplace date  if statefips == 36, c(l) lcolor(cranberry) ms(none) lwidth(medthick) lp(dash))
	,
	xtitle("") xmtick(##2)
	ytitle("Percent change in time spent in workplace")
	legend(order(22 "California" 23 "New York" 2 "Donor States") row(1) span) 
	yscale(r(-0.6 0.1))
	ylab(-0.5(0.2)0.1)
	yline(0, lcolor(black))
	graphregion(color(white)) bgcolor(white)
;
#delimit cr 
graph export "$projdir/out/figures/figure5b.pdf",  replace

*-------------------------------------------------------------------------------
* Figure 6 (left panel): 
* Standard Labor Supply
*-------------------------------------------------------------------------------

* Setup
clear all
set obs 1001

* Initialize parameter values
qui gen K = 1.5
qui gen M = 1.8
qui gen m = -0.20
qui gen L = 0.001*(_n - 1)

* Initial w(Ld)
qui gen w_Ld = 1 - K*L
qui replace w_Ld = . if w_Ld < 0

* Initial w(Ls)
qui gen w_Ls = (M*L+m)*(L <= 1)
qui replace w_Ls = . if L > 1
qui gen y_Ls = w_Ls

* Competitive equlibrium values
qui gen p_prime = (M+K*m)/(M + K)
qui gen w_prime= p_prime
qui gen L_prime = (1-m)/(M + K)
qui replace w_prime = . if L > L_prime
local Lprime = L_prime
local wprime = w_prime

* Initial monopsonist equilibrium values
qui gen MC_L_M = 2*M*L+m-0.15*M
qui replace MC_L_M = . if MC_L_M > 1
qui gen L_M_prime = (1-m+0.15*M)/(2*M + K)
qui gen p_M_prime = 2*M*L_M_prime+m - 0.15*M
qui gen w_M_prime = M*L_M_prime+m
local L_Mprime = L_M_prime
local w_Mprime = w_M_prime
local p_Mprime = p_M_prime
qui replace w_M_prime = . if L > L_M_prime

qui gen b = (2*M-0.15*2*M^2-0.15*M*K-2*M*m)/(1-m-0.15*2*M-0.15*K)
qui gen a = -0.15*b+0.15*M+m
qui gen MC_L_M_2 = (a + b*L)*(L <= 1)
qui replace MC_L_M_2 = . if MC_L_M_2 > 0.8

* Minimum wages
gen wmin_1 = 0.35*`wprime'+0.65*`w_Mprime'
gen L_wmin_1 = (wmin_1-m)/M
local wmin1 = wmin_1
local Lwmin_1 = L_wmin_1
qui replace wmin_1 = . if L > L_wmin_1

gen wmin_2 = 0.6*`wprime'+0.4*`p_Mprime'
gen L_wmin_2 = (1-wmin_2)/K
local wmin2 = wmin_2
local Lwmin_2 = L_wmin_2
qui replace wmin_2 = . if L > L_wmin_2

gen wmin_3 = `p_Mprime'+0.06
gen L_wmin_3 = (1-wmin_3)/K
local wmin3 = wmin_3
local Lwmin_3 = L_wmin_3
qui replace wmin_3 = . if L > L_wmin_3

* Restrict values for the graph
drop if L > 0.7 | L < 0.15
qui replace w_Ls= . if L > 0.7 | w_Ls > 0.8
qui replace MC_L_M = . if L > 0.7 | MC_L_M > 0.8

* Generate coordinates for change
local L_change_1_a = `L_Mprime' + 0.003
local L_change_1_b = `Lwmin_1' - 0.003
local w_change_1_a = `w_Mprime' + 0.005
local w_change_1_b = `wmin1' - 0.005

#delimit ;
	twoway
		(line w_Ld L, lcolor(black) lp(solid)) // L_d 
		(line w_Ls L, lcolor(black) lp(solid)) // L_s
		(line MC_L_M L, lcolor(black) lp(solid)) // MC
		(line w_prime L, lcolor(gs12) lp(dash) lwidth(thin)) //w_prime
		(line w_M_prime L, lcolor(gs8) lp(dash) lwidth(thin)) //w_M_prime
		(line wmin_1 L, lcolor(black) lp(solid) lwidth(medthick)) // wmin_1
		(pci 0 `Lprime' `wprime' `Lprime', lcolor(gs12) lp(dash) lwidth(thin)) // L_prime
		(pci 0 `L_Mprime' `p_Mprime' `L_Mprime', lcolor(gs8) lp(dash) lwidth(thin)) // L_M_prime
		(pci 0 `Lwmin_1' `wmin1' `Lwmin_1', lcolor(black) lp(dash) lwidth(thin)) // L_min_1
		(pcarrowi 0 0.15 0 0.7, lcolor(black) mcolor(black) lp(solid) msize(small)) // x-axis
		(pcarrowi 0 0.15 0.8 0.15, lcolor(black) mcolor(black) lp(solid) msize(small)) // y-axis
		(pcarrowi 0.05 `L_change_1_a' 0.05 `L_change_1_b', lcolor(black) mcolor(black) msize(small)) // chnage in L
		(pcarrowi `w_change_1_a' 0.18 `w_change_1_b' 0.18, lcolor(black) mcolor(black) msize(small)), // change in w
		text(0.05 0.66 "L{superscript:D}")
		text(0.72 0.54 "L{superscript:S}")
		text(0.72 0.36 "MC")
		text(0.8 0.136 "w", place(n))
		text(`wprime' 0.138 "w*")
		text(`w_Mprime' 0.136 "w")
		text(`w_Mprime' 0.142 "M", size(tiny) placement(ne))
		text(`wmin1' 0.136 "w")
		text(`wmin1' 0.142 "1", size(vsmall) placement(ne))
		text(`wmin1' 0.136 "__", size(vsmall) placement(s))
		text(-0.018 0.7 "L")
		text(-0.02 `Lprime' "L")
		text(-0.028 `Lprime' "*", placement(ne))
		text(-0.02 `L_Mprime' "L")
		text(-0.015 `L_Mprime' "M", size(tiny) placement(ne))
		text(-0.02 `Lwmin_1' "L")
		text(-0.02 `Lwmin_1' "1", size(vsmall) placement(ne))
		text(-0.02 `Lwmin_1' "__", size(vsmall) placement(s))
		graphregion(fcolor(white) margin(4 20 4 4) lcolor(white) ilcolor(white)) 
		plotregion(margin(0 0 0 0) lcolor(white))
		title("")
		xlabel(none)
		ylabel(none)
		legend(off);	
	#delimit cr

* Save	
graph export "$projdir/out/figures/figure6l.pdf", replace

*-------------------------------------------------------------------------------
* Figure 6 (left panel): 
* Flat Labor Supply
*-------------------------------------------------------------------------------

* Setup
clear all
set obs 1001

* Initialize parameter values
qui gen K = 1.5
qui gen M = 0.8
qui gen m = 0.07
qui gen L = 0.001*(_n - 1)

* Initial w(Ld)
qui gen w_Ld = 1 - K*L
qui replace w_Ld = . if w_Ld < 0

* Initial w(Ls)
qui gen w_Ls = (M*L+m)*(L <= 1)
qui replace w_Ls = . if L > 1
qui gen y_Ls = w_Ls

* Competitive equlibrium values
qui gen p_prime = (M+K*m)/(M + K)
qui gen w_prime= p_prime
qui gen L_prime = (1-m)/(M + K)
qui replace w_prime = . if L > L_prime
local Lprime = L_prime
local wprime = w_prime

* Initial monopsonist equilibrium values
qui gen MC_L_M = 2*w_Ls - m
qui replace MC_L_M = . if MC_L_M > 1
qui gen p_M_prime = (2*M+K*m)/(2*M + K)
qui gen L_M_prime = (1-m)/(2*M + K)
qui gen w_M_prime = M*L_M_prime+m
local L_Mprime = L_M_prime
local w_Mprime = w_M_prime
local p_Mprime = p_M_prime
qui replace w_M_prime = . if L > L_M_prime

qui gen c = 0
qui gen b = (2*M-0.15*2*M^2-0.15*M*K-2*M*m-2*M*c-c*K)/(1-m-0.15*2*M-0.15*K)
qui gen a = -0.15*b+0.15*M+m + c
qui gen MC_L_M_2 = (a + b*L)*(L <= 1)
qui replace MC_L_M_2 = . if MC_L_M_2 > 0.8

* Minimum wages
gen wmin_1 = 0.5*`wprime'+0.5*`w_Mprime'
gen L_wmin_1 = (wmin_1-m)/M
local wmin1 = wmin_1
local Lwmin_1 = L_wmin_1
qui replace wmin_1 = . if L > L_wmin_1

gen wmin_2 = 0.6*`w_Mprime' + 0.4*`wprime'
gen L_wmin_2 = (wmin_2-m)/M
local wmin2 = wmin_2
local Lwmin_2 = L_wmin_2
qui replace wmin_2 = . if L > L_wmin_2

gen wmin_3 = `p_Mprime'+0.15
gen L_wmin_3 = (1-wmin_3)/K
local wmin3 = wmin_3
local Lwmin_3 = L_wmin_3
qui replace wmin_3 = . if L > L_wmin_3

* Restrict values for the graph
drop if L > 0.7 | L < 0.15
qui replace w_Ls= . if L > 0.7 | w_Ls > 0.8
qui replace MC_L_M = . if L > 0.7 | MC_L_M > 0.8

* Generate coordinates for change
local L_change_1_a = `L_Mprime' + 0.003
local L_change_1_b = `Lwmin_2' - 0.003
local w_change_1_a = `w_Mprime' + 0.005
local w_change_1_b = `wmin2' - 0.005



#delimit ;
	twoway
		(line w_Ld L, lcolor(black) lp(solid)) // L_d 
		(line w_Ls L, lcolor(black) lp(solid)) // L_s
		(line MC_L_M_2 L, lcolor(black) lp(solid)) // MC_2
		(line w_prime L, lcolor(gs12) lp(dash) lwidth(thin)) //w_prime
		(line w_M_prime L, lcolor(gs8) lp(dash) lwidth(thin)) //w_M_prime
		(line wmin_2 L, lcolor(black) lp(solid) lwidth(medthick)) // wmin_2
		(pci 0 `Lprime' `wprime' `Lprime', lcolor(gs12) lp(dash) lwidth(thin)) // L_prime
		(pci 0 `L_Mprime' `p_Mprime' `L_Mprime', lcolor(gs8) lp(dash) lwidth(thin)) // L_M_prime
		(pci 0 `Lwmin_2' `wmin2' `Lwmin_2', lcolor(black) lp(dash) lwidth(thin)) // L_min_2
		(pcarrowi 0 0.15 0 0.7, lcolor(black) mcolor(black) lp(solid) msize(small)) // x-axis
		(pcarrowi 0 0.15 0.8 0.15, lcolor(black) mcolor(black) lp(solid) msize(small)) // y-axis
		(pcarrowi 0.05 `L_change_1_a' 0.05 `L_change_1_b', lcolor(black) mcolor(black) msize(small)) // change in L
		(pcarrowi `w_change_1_a' 0.18 `w_change_1_b' 0.18, lcolor(black) mcolor(black) msize(small)), // change in w
		text(0.05 0.66 "L{superscript:D}")
		text(0.54 0.66 "L{superscript:S}")
		text(0.74 0.41 "MC")
		text(0.8 0.136 "w", place(n))
		text(`wprime' 0.138 "w*")
		text(`w_Mprime' 0.136 "w")	
		text(`w_Mprime' 0.142 "M", size(tiny) placement(ne))
		text(`wmin2' 0.136 "w")
		text(`wmin2' 0.142 "1", size(vsmall) placement(ne))
		text(`wmin2' 0.136 "__", size(vsmall) placement(s))
		text(-0.018 0.7 "L")
		text(-0.02 `Lprime' "L")
		text(-0.028 `Lprime' "*", placement(ne))
		text(-0.02 `L_Mprime' "L")
		text(-0.015 `L_Mprime' "M", size(tiny) placement(ne))
		text(-0.02 `Lwmin_2' "L")
		text(-0.02 `Lwmin_2' "1", size(vsmall) placement(ne))
		text(-0.02 `Lwmin_2' "__", size(vsmall) placement(s))
		graphregion(fcolor(white) margin(4 20 4 4) lcolor(white)) 
		plotregion(margin(0 0 0 0) lcolor(white))
		title("")
		xlabel(none)
		ylabel(none)
		legend(off);	
	#delimit cr

* Save
graph export "$projdir/out/figures/figure6r.pdf", replace

* Close the log
log close sublog
